//
// File: xhseqr.cpp
//
// MATLAB Coder version            : 5.1
// C/C++ source code generated on  : 22-Feb-2021 03:11:30
//

// Include Files
#include "xhseqr.h"
#include "rt_nonfinite.h"
#include "xdlanv2.h"
#include "xzlarfg.h"
#include "rt_nonfinite.h"
#include <cmath>

// Function Definitions
//
// Arguments    : double h[16]
//                double z[16]
// Return Type  : int
//
namespace coder
{
  namespace internal
  {
    namespace lapack
    {
      int xhseqr(double h[16], double z[16])
      {
        double v[3];
        double aa;
        double ab;
        double ba;
        double bb;
        double rt1r;
        double rt2r;
        double s;
        double s_tmp;
        double tst;
        int hoffset;
        int i;
        int info;
        int iy;
        boolean_T exitg1;
        info = 0;
        v[0] = 0.0;
        v[1] = 0.0;
        v[2] = 0.0;
        h[2] = 0.0;
        h[3] = 0.0;
        h[7] = 0.0;
        i = 3;
        exitg1 = false;
        while ((!exitg1) && (i + 1 >= 1)) {
          int L;
          int b_i;
          int b_k;
          int its;
          int k;
          int m;
          int nr;
          boolean_T exitg2;
          boolean_T goto150;
          L = 1;
          goto150 = false;
          its = 0;
          exitg2 = false;
          while ((!exitg2) && (its < 301)) {
            boolean_T exitg3;
            k = i;
            exitg3 = false;
            while ((!exitg3) && (k + 1 > L)) {
              b_i = k + ((k - 1) << 2);
              ba = std::abs(h[b_i]);
              if (ba <= 4.0083367200179456E-292) {
                exitg3 = true;
              } else {
                nr = k + (k << 2);
                bb = std::abs(h[nr]);
                tst = std::abs(h[b_i - 1]) + bb;
                if (tst == 0.0) {
                  if (k - 1 >= 1) {
                    tst = std::abs(h[(k + ((k - 2) << 2)) - 1]);
                  }

                  if (k + 2 <= 4) {
                    tst += std::abs(h[nr + 1]);
                  }
                }

                if (ba <= 2.2204460492503131E-16 * tst) {
                  tst = std::abs(h[nr - 1]);
                  if (ba > tst) {
                    ab = ba;
                    ba = tst;
                  } else {
                    ab = tst;
                  }

                  tst = std::abs(h[b_i - 1] - h[nr]);
                  if (bb > tst) {
                    aa = bb;
                    bb = tst;
                  } else {
                    aa = tst;
                  }

                  s = aa + ab;
                  tst = 2.2204460492503131E-16 * (bb * (aa / s));
                  if ((4.0083367200179456E-292 > tst) || rtIsNaN(tst)) {
                    tst = 4.0083367200179456E-292;
                  }

                  if (ba * (ab / s) <= tst) {
                    exitg3 = true;
                  } else {
                    k--;
                  }
                } else {
                  k--;
                }
              }
            }

            L = k + 1;
            if (k + 1 > 1) {
              h[k + ((k - 1) << 2)] = 0.0;
            }

            if (k + 1 >= i) {
              goto150 = true;
              exitg2 = true;
            } else {
              if (its == 10) {
                nr = k + (k << 2);
                s = std::abs(h[nr + 1]) + std::abs(h[(k + ((k + 1) << 2)) + 2]);
                tst = 0.75 * s + h[nr];
                ab = -0.4375 * s;
                aa = s;
                ba = tst;
              } else if (its == 20) {
                s = std::abs(h[i + ((i - 1) << 2)]) + std::abs(h[(i + ((i - 2) <<
                  2)) - 1]);
                tst = 0.75 * s + h[i + (i << 2)];
                ab = -0.4375 * s;
                aa = s;
                ba = tst;
              } else {
                nr = i + ((i - 1) << 2);
                tst = h[nr - 1];
                aa = h[nr];
                ab = h[(i + (i << 2)) - 1];
                ba = h[i + (i << 2)];
              }

              s = ((std::abs(tst) + std::abs(ab)) + std::abs(aa)) + std::abs(ba);
              if (s == 0.0) {
                rt1r = 0.0;
                bb = 0.0;
                rt2r = 0.0;
                ba = 0.0;
              } else {
                tst /= s;
                aa /= s;
                ab /= s;
                ba /= s;
                bb = (tst + ba) / 2.0;
                tst = (tst - bb) * (ba - bb) - ab * aa;
                aa = std::sqrt(std::abs(tst));
                if (tst >= 0.0) {
                  rt1r = bb * s;
                  rt2r = rt1r;
                  bb = aa * s;
                  ba = -bb;
                } else {
                  rt1r = bb + aa;
                  rt2r = bb - aa;
                  if (std::abs(rt1r - ba) <= std::abs(rt2r - ba)) {
                    rt1r *= s;
                    rt2r = rt1r;
                  } else {
                    rt2r *= s;
                    rt1r = rt2r;
                  }

                  bb = 0.0;
                  ba = 0.0;
                }
              }

              m = i - 1;
              exitg3 = false;
              while ((!exitg3) && (m >= k + 1)) {
                nr = m + ((m - 1) << 2);
                tst = h[nr];
                s_tmp = h[nr - 1];
                aa = s_tmp - rt2r;
                s = (std::abs(aa) + std::abs(ba)) + std::abs(tst);
                ab = tst / s;
                nr = m + (m << 2);
                v[0] = (ab * h[nr - 1] + (s_tmp - rt1r) * (aa / s)) - bb * (ba /
                  s);
                tst = h[nr];
                v[1] = ab * (((s_tmp + tst) - rt1r) - rt2r);
                v[2] = ab * h[nr + 1];
                s = (std::abs(v[0]) + std::abs(v[1])) + std::abs(v[2]);
                v[0] /= s;
                v[1] /= s;
                v[2] /= s;
                if (m == k + 1) {
                  exitg3 = true;
                } else {
                  b_i = m + ((m - 2) << 2);
                  if (std::abs(h[b_i - 1]) * (std::abs(v[1]) + std::abs(v[2])) <=
                      2.2204460492503131E-16 * std::abs(v[0]) * ((std::abs(h[b_i
                         - 2]) + std::abs(s_tmp)) + std::abs(tst))) {
                    exitg3 = true;
                  } else {
                    m--;
                  }
                }
              }

              for (b_k = m; b_k <= i; b_k++) {
                int j;
                nr = (i - b_k) + 2;
                if (3 < nr) {
                  nr = 3;
                }

                if (b_k > m) {
                  hoffset = (b_k + ((b_k - 2) << 2)) - 1;
                  for (j = 0; j < nr; j++) {
                    v[j] = h[j + hoffset];
                  }
                }

                tst = v[0];
                bb = reflapack::xzlarfg(nr, &tst, v);
                v[0] = tst;
                if (b_k > m) {
                  h[(b_k + ((b_k - 2) << 2)) - 1] = tst;
                  h[b_k + ((b_k - 2) << 2)] = 0.0;
                  if (b_k < i) {
                    h[b_k + 1] = 0.0;
                  }
                } else {
                  if (m > k + 1) {
                    h[b_k - 1] *= 1.0 - bb;
                  }
                }

                s = v[1];
                aa = bb * v[1];
                if (nr == 3) {
                  rt1r = v[2];
                  tst = bb * v[2];
                  for (j = b_k; j < 5; j++) {
                    iy = b_k + ((j - 1) << 2);
                    ab = (h[iy - 1] + s * h[iy]) + rt1r * h[iy + 1];
                    h[iy - 1] -= ab * bb;
                    h[iy] -= ab * aa;
                    h[iy + 1] -= ab * tst;
                  }

                  if (b_k + 3 < i + 1) {
                    b_i = b_k + 2;
                  } else {
                    b_i = i;
                  }

                  for (j = 0; j <= b_i; j++) {
                    iy = j + ((b_k - 1) << 2);
                    nr = j + (b_k << 2);
                    hoffset = j + ((b_k + 1) << 2);
                    ab = (h[iy] + s * h[nr]) + rt1r * h[hoffset];
                    h[iy] -= ab * bb;
                    h[nr] -= ab * aa;
                    h[hoffset] -= ab * tst;
                  }

                  for (j = 0; j < 4; j++) {
                    iy = j + ((b_k - 1) << 2);
                    nr = j + (b_k << 2);
                    hoffset = j + ((b_k + 1) << 2);
                    ab = (z[iy] + s * z[nr]) + rt1r * z[hoffset];
                    z[iy] -= ab * bb;
                    z[nr] -= ab * aa;
                    z[hoffset] -= ab * tst;
                  }
                } else {
                  if (nr == 2) {
                    for (j = b_k; j < 5; j++) {
                      iy = b_k + ((j - 1) << 2);
                      tst = h[iy - 1];
                      ab = tst + s * h[iy];
                      h[iy - 1] = tst - ab * bb;
                      h[iy] -= ab * aa;
                    }

                    for (j = 0; j <= i; j++) {
                      iy = j + ((b_k - 1) << 2);
                      nr = j + (b_k << 2);
                      ab = h[iy] + s * h[nr];
                      h[iy] -= ab * bb;
                      h[nr] -= ab * aa;
                    }

                    for (j = 0; j < 4; j++) {
                      iy = j + ((b_k - 1) << 2);
                      tst = z[iy];
                      nr = j + (b_k << 2);
                      ab = tst + s * z[nr];
                      z[iy] = tst - ab * bb;
                      z[nr] -= ab * aa;
                    }
                  }
                }
              }

              its++;
            }
          }

          if (!goto150) {
            info = i + 1;
            exitg1 = true;
          } else {
            if ((L != i + 1) && (L == i)) {
              m = i << 2;
              b_i = i + m;
              s = h[b_i - 1];
              b_k = (i - 1) << 2;
              nr = i + b_k;
              rt1r = h[nr];
              tst = h[b_i];
              reflapack::xdlanv2(&h[(i + ((i - 1) << 2)) - 1], &s, &rt1r, &tst,
                                 &aa, &ab, &bb, &ba, &s_tmp, &rt2r);
              h[b_i - 1] = s;
              h[nr] = rt1r;
              h[b_i] = tst;
              if (4 > i + 1) {
                nr = 2 - i;
                iy = i + ((i + 1) << 2);
                hoffset = iy - 1;
                for (k = 0; k <= nr; k++) {
                  tst = s_tmp * h[hoffset] + rt2r * h[iy];
                  h[iy] = s_tmp * h[iy] - rt2r * h[hoffset];
                  h[hoffset] = tst;
                  iy += 4;
                  hoffset += 4;
                }
              }

              if (i - 1 >= 1) {
                hoffset = b_k;
                iy = m;
                for (k = 0; k <= i - 2; k++) {
                  tst = s_tmp * h[hoffset] + rt2r * h[iy];
                  h[iy] = s_tmp * h[iy] - rt2r * h[hoffset];
                  h[hoffset] = tst;
                  iy++;
                  hoffset++;
                }
              }

              tst = s_tmp * z[b_k] + rt2r * z[m];
              z[m] = s_tmp * z[m] - rt2r * z[b_k];
              z[b_k] = tst;
              tst = z[m + 1];
              aa = z[b_k + 1];
              z[m + 1] = s_tmp * tst - rt2r * aa;
              z[b_k + 1] = s_tmp * aa + rt2r * tst;
              tst = z[m + 2];
              aa = z[b_k + 2];
              z[m + 2] = s_tmp * tst - rt2r * aa;
              z[b_k + 2] = s_tmp * aa + rt2r * tst;
              tst = z[m + 3];
              aa = z[b_k + 3];
              z[m + 3] = s_tmp * tst - rt2r * aa;
              z[b_k + 3] = s_tmp * aa + rt2r * tst;
            }

            i = L - 2;
          }
        }

        h[3] = 0.0;
        return info;
      }
    }
  }
}

//
// File trailer for xhseqr.cpp
//
// [EOF]
//
